/*
 *  lib-smv: Matrix - Vector Library targeted to Structural Mechanics Problems
 *  Copyright (C) 2006-2008  St.Brcic, Lj.Zugic
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
/** @file iArrays.cpp
 *  @brief Implementation of iVector and iMatrix classes of int numbers 
 */

#include "iArrays.h"

namespace smv_lib
{

const int iM_DEFAULT = 2;  // default size
//Defining ARRAYS_DEBUG turns on bounds checking and size checking
//#define ARRAYS_DEBUG

/*==========================================
      iVector class definitions
============================================*/
//! default constructor
iVector::iVector()                 // Vector of integers
{
  data = new int[iM_DEFAULT];
  assert(data != NULL);
  len = iM_DEFAULT;
  Name = "iVECTOR";
  isNameGiven = false;
}

//! constructor with given length and elements
iVector::iVector(int n, int d) : len(n)
{
  data = new int[len];
  assert(data != NULL);
  for (int i = 0; i < len; i++)
    data[i] = d;
  Name = "iVECTOR";
  isNameGiven = false;
}

//! constructor with a length and a name
iVector::iVector(int n, string name) : len(n), Name(name)
{
  data = new int[n];
  assert(data != NULL);
  int d = 0;
  for (int i = 0; i < n; i++)
    data[i] = d;
  if (Name == "") {
    Name = "iVECTOR";
    isNameGiven = false;
  }
  isNameGiven = true;
}

static char iAlphabet[] = "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ";
//! constructor from a file
//-------------------------------------------------------------------------------------
/**
** \brief Construct an iVector from a file
**
**    Format of the Vector file - it contains two parts: Header block and Data block.
**    HEADER BLOCK: 
**
**    (1) One or more optional COMMENT lines, where each comment line begins with '!'
**    
**    (2) Optional line with a NAME of the Vector. If given, it must be written as:
**    N= (string), 
**    where the string in the name must begin with a letter (a-w,A-W).
**    NOTE: The string must not contain space or tabs.
**    If the Name is not given, default name VECTOR is assigned.
**
**    (3) Number of elements 'n' and number of entries in a line 'l'
**    (for small vectors n and l are equal - all in one line).
**    Both entries 'n' and 'l' must be given!!
**    n= (int)  l= (int)
**           
**    The end of the Header block is given as 3 colons:
**    :::
**    DATA BLOCK:
**    numeric entris are given in a row 
**    (entries are separated by a blanko)
**
**    NOTE: For longer vectors (say, n > 15-20), multiline input is neccessary 
**    
**/
//-------------------------------------------------------------------------------------
iVector::iVector(const char* file)
{
  const int BSIZE = 120;
  FILE *from;
  int nele  = 0;  // total number of elements in a vector
  int mult  = 0;  // number of elements in a single line

  if ((from = fopen(file,"r")) == 0) {
    cerr << "\nCannot open " << file << endl;
    exit(1);
  }
  char buf[BSIZE], *cp, *cp2;
  int efound = 0, mfound = 0, nfound = 0, colonsfound = 0;

  //--------------------------------------------- read the HEADER BLOCK
  while( fgets(buf, BSIZE, from) ) {
    //--- read the comment lines
    if( ( cp = strpbrk(buf,"!")) != NULL )
      *cp = '\0';
    //---- read the Name of a Vector
    if(( cp = strpbrk(buf,"N")) != NULL )
      if (( cp2 = strpbrk(cp, "=")) != NULL )
        if (( cp = strpbrk( cp2, iAlphabet)) != NULL ) {
          Name = cp;
          isNameGiven = true;
          nfound++;
        }
    //--- read the lenth of a Vector
    if(( cp = strpbrk(buf,"n")) != NULL )
      if (( cp2 = strpbrk(cp, "=")) != NULL )
        if (( cp = strpbrk( cp2, "0123456789")) != NULL ) {
          nele = atoi(cp);
          efound++;
        }
    //--- read the number of entries in a line
    if(( cp = strpbrk(buf,"l")) != NULL )
      if (( cp2 = strpbrk(cp, "=")) != NULL )
        if (( cp = strpbrk( cp2, "0123456789")) != NULL ) {
          mult = atoi(cp);
          mfound++;
        }
    //--- termination of the Header block
    if ( strstr(buf,":::") != NULL ) {
      colonsfound++;
      break; //... out of "while" loop
    }
  }
  if ( !efound || ! mfound  || !colonsfound ) {
    cerr << "\nThe input file is not OK \n" << endl;
    exit(1);
  }
  if (!nfound) {               // no Vector name is given
    Name = "iVECTOR";         // assign the default name
    isNameGiven = false;
  }
  //------------------------------------------ make space for the vector
  data = new int[nele];
  assert (data != NULL);
  len  = nele;
  int i = 0, num;
  // number of line to read
  int No_lines, resto;
  div_t Res;
  Res.quot = 1;
  Res.rem = 0;
  if (nele > mult)
    Res = div(nele, mult);
  No_lines = Res.quot;
  resto    = Res.rem;
  //-------------------------------------------- read the DATA BLOCK
  for (int j = 0; j < No_lines; j++) {
    for (int k = 0; k < mult; k++) {
      char nb[20];
      fscanf(from, "%s", nb);         // scan for space-delimited string
      num = atoi(nb);                 // convert string to an int value
      i = j * mult + k;               // position in the vector
      setData(i,num);                 // assign data to a vector
      if (i >= len)
        break;
    }
  }
  // close the file
  fclose(from);
  if (Name != "iVECTOR") {
    isNameGiven = true;
  }
}

//! copy constructor
iVector::iVector(const iVector& original)
{
  data = new int[len = original.getLen()];
  assert (data != NULL);
  for(int i=0; i<len; i++)
    data[i] = original.data[i];
  Name = original.getName();
  if (Name == "iVECTOR")
    isNameGiven = false;
  else
    isNameGiven = true;
}

//! destructor
iVector::~iVector()
{
  delete [] data;
}

//--- operator functions
//! assignment operator
iVector& iVector::operator = (const iVector &original)
{
  if(this != &original) {
    delete [] data;
    data= new int[len = original.getLen()];
    assert (data != NULL);
    for(int i=0; i<len; i++)
      data[i] = original.data[i];
    Name = original.getName();
    if (Name == "iVECTOR")
      isNameGiven = false;
    else
      isNameGiven = true;
  }
  return *this;
}

//! operator []: set an element to the iVector
int& iVector::operator [] (int i)
{
  assert(i>=0 && i <len);

  return data[i];
}

//! operator []: get an element from the iVector
const int& iVector::operator [] (int i) const
{
  assert(i>=0 && i <len);

  return data[i];
}

const int iVector::getData(int i) const
{
  assert(i>=0 && i <len);

  return data[i];
}

void iVector::setData(int i, int d)
{
  assert(i>=0 && i <len);

  data[i] = d;
}

void iVector::setName(string name)
{
  Name = name;
  isNameGiven = true;
}

//! iVector b = self + v
iVector iVector::operator + (const iVector& v)
{
  assert(len == v.getLen());

  iVector sum(len);
  for(int i=0; i<len; i++)
    sum.data[i] = data[i] + v.data[i];
  return (sum);
}

//! iVector b = self - v
iVector iVector::operator - (const iVector& v)
{
  assert(len == v.getLen());

  iVector sum(len);
  for(int i=0; i<len; i++)
    sum.data[i] = data[i] - v.data[i];
  return (sum);
}

//! Operator '*' with an int: iVector v = self * s;
iVector iVector::operator * ( int s )
{
  iVector ans( len );
  for ( int i = 0; i < len; i++ )
    ans.data[ i ] = data[ i ] * s ;
  return ( ans );
}

//! scalar product with a vector: s = self * v
const int iVector::operator * ( const iVector& v ) const
{
  assert( len == v.getLen() );

  int s = 0;
  for ( int i = 0; i < len; i++ )
    s += data[i] * v.data[ i ];
  return s;
}

//! Operator '+=' iVector self = self + v
iVector& iVector::operator += ( const iVector& v )
{
  assert( len == v.getLen() );

  for ( int i = 0; i < len; i++ )
    data[ i ] += v.data[ i ];
  return ( *this );
}

//! Operator '-=' iVector self = self - v
iVector& iVector::operator -= ( const iVector& v )
{
  assert( len == v.getLen() );

  for ( int i = 0; i < len; i++ )
    data[ i ] -= v.data[ i ];
  return ( *this );
}

//! Operator '*=' iVector self = self * s
iVector& iVector::operator *= ( int s )
{
  for ( int i = 0; i < len; i++ )
    data[ i ] *= s;
  return ( *this );
}

//------------------------- umesto njih mogu da se koriste Add i Sub
//! iVector b = self + v
iVector iVector::Add(const iVector& v)
{
  assert(len == v.getLen());

  iVector sum(len);
  for(int i=0; i<len; i++)
    sum.data[i] = data[i] + v.data[i];
  return (sum);
}

//! iVector b = self - v
iVector iVector::Sub(const iVector& v)
{
  assert(len == v.getLen());

  iVector sum(len);
  for(int i=0; i<len; i++)
    sum.data[i] = data[i] - v.data[i];
  return (sum);
}

//! Operator '*' with an int: iVector v = self * s;
iVector iVector::Dot( int s )
{
  iVector ans( len );
  for ( int i = 0; i < len; i++ )
    ans.data[ i ] = data[ i ] * s;
  return ( ans );
}

//! scalar product with a vector: s = self * v
const int iVector::Dot( const iVector& v ) const
{
  assert( len == v.getLen() );

  int s = 0;
  for ( int i = 0; i < len; i++ )
    s += data[i] * v.data[ i ];
  return s;
}

//! outer product of self and vector v
iMatrix iVector::OuterDot( const iVector& v )
{
  int m = len;
  int n = v.getLen();
  iMatrix ans( m, n );
  int A_ij;
  for ( int i = 0; i < m; i++ ) {
    for ( int j = 0; j < n; j++ ) {
      A_ij = data[ i ] * v.data[ j ];
      ans.setElement( i, j, A_ij );
    }
  }
  return ( ans );
}

//! true if self == v
bool iVector::operator == (const iVector& v)
{
  if(len != v.getLen())
    return false;
  for(int i=0; i<len; i++)
    if(data[i]!=v.data[i])
      return false;
  if (Name != v.getName())       // names are not so important
    cerr << "\nVector names are NOT equal!" << endl;
  //vectors are equal
  return true;
}

//! Add sub-vector into the vector (add to the old content)
iVector& iVector::addSubVector( int i, const iVector& Sub )
{
  int r = Sub.getLen();
  // check if Sub-vector fits into this vector
  assert( (i+r) <= len );
  int ai;
  for ( int ii = 0; ii < r; ii++ ) {
    ai = this->getData( i + ii );
    this->setData( i + ii, ai + Sub.getData( ii ) );
  }
  return *this;
}

//! Take sub-vector from a given vector
iVector iVector::getSubVector( int i, int subLen )
{
  iVector result( subLen );
  // check if subvector belongs to this vector
  assert( (i+subLen) <= len );
  int ai;
  for ( int ii = 0; ii < subLen; ii++ ) {
    ai = this->getData( i + ii );
    result.setData( ii, ai );
  }
  return result;
}

// ------------------- arithmetic operations
//! change the size of a Vector
void iVector::Resize(int n)
{
  delete [] data;
  data = new int[len = n];
  assert(data != NULL);
}

//! make all elements equal to zero
void iVector::Zero(void)
{
  for(int i=0; i<len; i++)
    data[i]=0;
}

//! max element in a Vector
const int iVector::Max(void) const
{
  int max = data[0];
  for(int i=1; i<len; i++)
    if(data[i] > max)
      max = data[i];
  return max;
}

//! min element in a Vector
const int iVector::Min(void) const
{
  int min = data[0];
  for(int i=1; i<len; i++)
    if(data[i] < min)
      min = data[i];
  return min;
}

//! mean value of all elements
const double iVector::Mean(void) const
{
  double sum = 0.0;
  for(int i=0; i<len; i++)
    sum += data[i];
  return sum / len;
}

//! norm of Self vector (intensity)
const double iVector::Norm(void) const
{
  int s = 0;
  for(int i=0; i<len; i++)
    s+= data[i] * data[i];
  return sqrt((double)s);
}

//----------------- console I/O
//! print vector on screen in a single line
void  iVector::rPrint(void)  //print entries on a single line
{
  cout << "Vector: " << Name << " = ";
  cout << "(";
  for(int i=0; i<len-1; i++)
    cout << data[i] << ",";
  cout << data[len-1] << ")" << endl;
  return;
}

//! print a vector on screen
void iVector::Print(void)
{
  cout << "Vector:  " << Name << endl;
  cout << "Length = " << len  << endl;
  for(int i=0; i<len; i++)
    cout << data[i] << endl;
  return;
}

//------------------- File I/O
//
//! read vector from predefined file
iVector iVector::Read(void)
{
  string FileName;
  FileName = Name + ".txt";
  iVector ans(FileName.c_str());
  return (ans);
}

//! read form a given file
iVector iVector::ReadFrom(const char* file)
{
  iVector ans(file);
  return (ans);
}

//! write to a predefined file
void iVector::Write()
{
  int L = 20;   // default number of entries in a line
  if (len < L)
    L = len;
  string FileName;
  FileName = Name + ".txt";
  this->WriteTo(FileName.c_str(),L);
}

//! write to a given file
void iVector::WriteTo(const char* file, int L)
{
  int No_Lines;          // number of lines for ONE ROW
  int resto;             // left-over in the last line for the row
  div_t Res;
  Res.quot = 1;
  Res.rem = 0;
  if (len <= L)
    L = len;
  else if (len > L)                    // if number of cols is > number of entries in a row
    Res = div(len, L);
  No_Lines = Res.quot;                 // No_Lines with full number of mult entries
  resto    = Res.rem;                  // ONE line with resto entries
  // all for just ONE ROW of data
  int J;
  string VectorName = Name;
  //-- output to file
  out.open(file,ios::out);
  if (!out)
    exit(2);
  //--- write the HEADER BLOCK
  out << "!================== iVector(" << len << ") ====================" << endl;
  out << "N= " << VectorName << endl;
  out << "n= " << len <<  "  l= " << L  << endl;
  out << ":::" << endl;
  //--- write the DATA BLOCK
  int kmax;
  for (int j = 0; j <= No_Lines; j++) {
    kmax = L;
    if (j == No_Lines)
      kmax = resto;
    for (int k = 0; k < kmax; k++) {
      J = j * L + k;
      out << this->getData(J) << "   ";
    }
    out << endl;
  }
  out.close();
}

//! global iVector operator functions
//=======================================================================
//! Addition of two vectors: c = a + b
iVector operator + (const iVector& a, const iVector& b)
{
  // check compatibility of vectors
  int n = a.getLen();
  assert(n == b.getLen());
  iVector c(n);
  int ci;
  for (int i=0; i<n; i++) {
    ci = a.getData(i) + b.getData(i);
    c.setData(i,ci);
  }
  return (c);
}

//! Subtraction of two vectors: c = a - b
iVector operator - (const iVector& a, const iVector& b)
{
  // check compatibility of vectors
  int n = a.getLen();
  assert(n == b.getLen());
  iVector c(n);
  int ci;
  for (int i=0; i<n; i++) {
    ci = a.getData(i) - b.getData(i);
    c.setData(i,ci);
  }
  return (c);
}

//! dot product s = v1 * v2
int Dot(const iVector& v1, const iVector& v2)
{
  // check vector compatibility
  assert(v1.getLen() == v2.getLen());
  int s = 0;
  for(int i=0; i<v1.getLen(); i++)
    s += v1.getData(i) * v2.getData(i);
  return s;
}

//! dot product of two vectors: d = v1 * v2
int operator * (const iVector& v1, const iVector& v2)
{
  // check vector compatibility
  assert(v1.getLen() == v2.getLen());
  int s = 0;
  for(int i=0; i<v1.getLen(); i++)
    s += v1.getData(i) * v2.getData(i);
  return s;
}

//! Product with int scalar: Vector b = s * v
iVector operator * (int x, const iVector& v)
{
  iVector ans(v.getLen());
  for(int i=0; i<v.getLen(); i++)
    ans.setData(i, x*v[i]);
  return (ans);
}

//! Product with int scalar: Vector b = v * s
iVector operator * (const iVector& v, int x)
{
  iVector ans(v.getLen());
  for(int i=0;i<v.getLen();i++)
    ans.setData(i, x*v[i]);
  return (ans);
}

//! norm of a vector v
double Norm(const iVector& v)
{
  int s = 0;
  for(int i=0; i< v.getLen(); i++)
    s+= v.getData(i) * v.getData(i);
  return sqrt((double)s);
}

//! output stream
ostream& operator << (ostream& s, const iVector& v)
{
  s << "(";
  for(int i=0; i<v.getLen()-1; i++)
    s << v.getData(i) << ", ";
  s << v.getData(v.getLen()-1) <<  ")"  <<endl;
  return s;
}

//================================================== end of iVECTOR CLASS
/*===============================================
          iMatrix class definitions
===============================================*/
//
//----------------- constructors
//--- empty constructor - makse a matrix 2 x 2 with all zeroes
//    iM_DEFAULT = 2 (can be changed)
//
//! default constructor
iMatrix::iMatrix()
{
  m = new iVector[iM_DEFAULT];
  assert(m != NULL);
  rows = cols = iM_DEFAULT;
  iVector v(cols);
  for(int i=0; i<rows; i++)
    m[i]= v;
  Name = "iMATRIX";
  isNameGiven = false;
}

//--- iMatrix r x c, all elements are equal to dat
//
//! constructor with given dimensions and equal data
iMatrix::iMatrix(int r, int c, int dat) : rows(r), cols(c)
{
  m = new iVector[rows];
  assert(m != NULL);
  iVector v(cols,dat);
  for(int i=0; i<r; i++)
    m[i] = v;
  Name = "iMATRIX";
  isNameGiven = false;
}

//! constructor with given dimensions and a name
iMatrix::iMatrix(int r, int c, string name) : rows(r), cols(c), Name(name)
{
  m = new iVector[r];
  assert(m != NULL);
  int d = 0;
  iVector v(cols,d);
  for(int i=0; i<r; i++)         // Matrix (r x c) with all elements set to zero
    m[i] = v;
  if (Name == "") {
    Name = "iMATRIX";
    isNameGiven = false;
  }
  isNameGiven = true;
}

//--------------------------------------------------------------------------------------------
/**
** \brief Construct an iMatrix from a file
**
**    Format of the matrix file - it contains two parts: Header block and Data block
**    HEADER BLOCK:
**
**    (1) One or more optional COMMENT lines, where each comment line begins with '!'
**    
**    (2) Optional line with a NAME of the Matrix. If given, it must be written as:
**    N= (string) 
**
**    where the string in the name must begin with a letter (a-w,A-W).
**    NOTE: The string must not contain space or tabs.
**    If the Name is not given, default name MATRIX is assigned.
**
**    (3) Number of ROWS 'r', number of COLUMNS 'c' and number of elements in one line 'l'
**    r= (int)  c= (int)  l= (int)
**
**    (for small number of columns, c = l, i.e. each row is in one line)
**    NOTE - all three entries, r, c, and l must be given!!
**
**    The end of the Header part is given as 3 colons:
**    :::
**
**    DATA BLOCK:
**    numeric entris are given row by row (number of entries in a row is <= l)
**    (entries are separated by a blanko)
**
**    NOTE: For larger matrices (say, cols > 15-20), multiline input for each row is neccessary 
**    
**/
//----------------------------------------------------------------------------------------------
//! constructor from a file
iMatrix::iMatrix(const char* file)
{
  const int BSIZE = 120;
  FILE *from;
  int row = 0, col = 0;       // number of rows and columns
  int mul = 0;                // number of entries in one line

  if ((from = fopen(file,"r")) == 0) {
    cerr << "\nCannot open " << file << endl;
    exit(1);
  }
  char buf[BSIZE], *cp, *cp2;
  int rfound = 0, cfound = 0, mfound = 0, nfound = 0, colonsfound = 0;  // control flags

  //---------------------------------------------------- read the HEADER BLOCK
  while( fgets(buf, BSIZE, from) ) {
    //--- read the comment line
    if( ( cp = strpbrk(buf,"!")) != NULL )
      *cp = '\0';
    //---- read the Name of a Vector
    if(( cp = strpbrk(buf,"N")) != NULL )
      if (( cp2 = strpbrk(cp, "=")) != NULL )
        if (( cp = strpbrk( cp2, iAlphabet)) != NULL ) {
          Name = cp;
          isNameGiven = true;
          nfound++;
        }
    //--- read the number of rows
    if(( cp = strpbrk(buf,"r")) != NULL )
      if (( cp2 = strpbrk(cp, "=")) != NULL )
        if (( cp = strpbrk( cp2, "0123456789")) != NULL ) {
          row = atoi(cp);
          rfound++;
        }
    //--- read number of columns
    if( ( cp = strpbrk(buf,"c") ) != NULL )
      if ( ( cp2 = strpbrk(cp, "=")) != NULL )
        if ( ( cp = strpbrk(cp2, "0123456789")) != NULL ) {
          col = atoi(cp);
          cfound++;
        }
    //--- read number of entries in a line
    if( ( cp = strpbrk(buf,"l") ) != NULL )
      if ( ( cp2 = strpbrk(cp, "=")) != NULL )
        if ( ( cp = strpbrk(cp2, "0123456789")) != NULL ) {
          mul = atoi(cp);
          mfound++;
        }
    //--- termination of the Header block
    if ( strstr(buf,":::") != NULL ) {
      colonsfound++;
      break; //... out of "while" loop
    }
  }
  if ( !rfound || !cfound || !mfound || !colonsfound ) { // if the header block is finished
    cerr << "\nThe input file is NOT OK \n" << endl;
    exit(1);
  }
  if (!nfound) {              // no Matrix name is given
    Name = "iMATRIX";         // assign the default name
    isNameGiven = false;
  }
  //------------------------------------------------- make space for the matrix
  m = new iVector[row];             // array of row vectors
  assert(m != NULL);
  iVector* v = new iVector(col);    // each vector with col elements
  assert(v != NULL);
  rows = row;
  cols = col;
  int I = 0, data;
  // number of input lines for each matrix row (No_Lines) and
  // number of remaining elements in the last line for a row (resto)
  // --- for one row of the matrix
  int No_Lines, resto;
  div_t Res;
  Res.quot = 1;
  Res.rem = 0;
  if (col > mul)                   // if number of cols is > number of entries in a row
    Res = div(col, mul);
  No_Lines = Res.quot;             // No_Lines with full number of mult entries
  resto    = Res.rem;              // ONE line with resto entries, all for just ONE ROW of data
  //-------------------------------------- read the DATA BLOCK
  int kmax;
  for (int i = 0; i < rows; i++) {
    for (int j = 0; j <= No_Lines; j++) {
      kmax = mul;
      if (j == No_Lines)                 // proper number of data
        kmax = resto;                    // for the last line in a row
      for (int k = 0; k < kmax; k++) {     // read a row from a file
        char nb[20];
        fscanf(from, "%s", nb);          // scan for space-delimited string
        data = atoi(nb);                 // convert string to an int value
        I = j * mul + k;
        v->setData(I,data);              // assign data to each row
      }
    }  // end of data for one row
    m[i] = *v;
  }
  fclose(from);
  if (Name != "iMATRIX") {
    //Name.TrimRight();
    isNameGiven = true;
  }
}

//! copy constructor
iMatrix::iMatrix(const iMatrix& s)
{
  rows = s.rows;
  m = new iVector[rows];
  assert(m != NULL);
  cols = s.cols;
  for(int i=0; i<rows; i++) {
    m[i] = s.m[i];
  }
  Name = s.getName();
  if (Name == "iMATRIX")
    isNameGiven = false;
  else
    isNameGiven = true;
}

//! destructor
iMatrix::~iMatrix()
{
  delete [] m;
}

//! assignment operator
iMatrix& iMatrix::operator = (const iMatrix& s)
{
  if(this != &s) {
    delete []m;
    rows = s.getRows();
    cols = s.getCols();
    m = new iVector[rows];
    assert(m != NULL);
    for(int i=0; i<rows; i++)
      this->m[i] = s.getRowData(i);
    Name = s.getName();
    if (Name == "iMATRIX")
      isNameGiven = false;
    else
      isNameGiven = true;
  }
  return *this;
}

//! get element (i,j) - zero based
const int iMatrix::getElement(int i, int j) const
{
  assert(i>=0 && i<rows);
  assert(j>=0 && j<cols);

  int val = this->m[i].getData(j);
  return (val);
}

//! set the value 'val' to row 'i' and column 'j'
void iMatrix::setElement(int i, int j, int val)
{
  assert(i>=0 && i<rows);
  assert(j>=0 && j<cols);

  this->m[i].setData(j,val);
}

//! set a name to the Matrix
void iMatrix::setName(string name)
{
  Name = name;
  isNameGiven = true;
}

//! get elements in a row 'i' of a Matrix
const iVector iMatrix::getRowData(int i) const
{
  assert(i>=0 && i<rows);

  return m[i];
}

//! get elements in a column 'j' of a Matrix
const iVector iMatrix::getColData(int j) const
{
  assert(j>=0 && j<cols);

  int i, Rows = this->getRows();
  int val;
  iVector v(Rows);
  for (i = 0; i < Rows; i++) {
    val = m[i].getData(j);
    v.setData(i,val);
  }
  return (v);
}

//! set all data in a column 'j' with Vector 'v'
void iMatrix::setColData(int j, const iVector& v)
{
  assert(j>0 && j<cols);
  int Vrows = v.getLen();
  assert(rows == Vrows);

  int i, Rows = this->getRows();
  int val;
  for (i = 0; i < Rows; i++) {
    val = v.getData(i);
    m[i].setData(j,val);
  }
}

//! set all data in a row 'i' with Vector 'v'
void iMatrix::setRowData(int i, const iVector& v)
{
  assert(i>=0 && i<rows);
  int Cols = v.getLen();
  assert(cols == Cols);

  this->m[i] = v;
}

//------------------- arithmetic operations
//! Sum of elements in a row 'i'
//
const int iMatrix::RowSum(int i) const
{
  assert(i>=0 && i<rows);

  int sum = 0;
  for(int k=0; k<cols; k++)
    sum += m[i][k];
  return sum;
}

//! Sum of elements in a column 'j'
//
const int iMatrix::ColSum(int j) const
{
  assert(j>=0 && j<cols);

  int sum = 0;
  for(int k=0; k<rows; k++)
    sum += m[k][j];
  return sum;
}

//! Matrix is assigned all zero elements
//
void iMatrix::Zero(void)
{
  for(int i=0; i<rows; i++) 
    for(int j=0;j<cols;j++)
      m[i][j] =0 ;
}

//! Matrix is transformed into an identity matrix
//
void iMatrix::Identity(void)
{
  assert(rows == cols);

  Zero();
  for(int i=0; i<rows; i++)
    m[i][i]=1;
}

//! Matrix is transformed into a diagonal matrix (all elements equal to 'd')
//
void iMatrix::Diagonal(int d)
{
  assert(rows == cols);

  Zero();
  for(int i=0; i<rows; i++)
    m[i][i] = d;
}

//! Matrix is transformed into a diagonal matrix (elements equal to iVector 'v')
//
void iMatrix::Diagonal(const iVector& v)
{
  assert(rows == cols);

  int d_i;
  Zero();
  for(int i=0; i<rows; i++) {
    d_i = v.getData(i);
    m[i][i] = d_i;
  }
}

//! Insert sub-matrix into the matrix (replace the old content)
iMatrix& iMatrix::insertSubMatrix(int i, int j, const iMatrix& Sub)
{
  int r = Sub.getRows();
  int c = Sub.getCols();
  // check if Sub-matrix fits into this matrix
  assert( (i+r) <= rows || (j+c) <= cols );
  for (int ii=0; ii<r; ii++) {
    for (int jj=0; jj<c; jj++)
      this->setElement(i+ii,j+jj, Sub.getElement(ii,jj));
  }
  return *this;
}

//! Add sub-matrix into the matrix (add to the old content)
iMatrix& iMatrix::addSubMatrix(int i, int j, const iMatrix& Sub)
{
  int r = Sub.getRows();
  int c = Sub.getCols();
  // check if Sub-matrix fits into this matrix
  assert( (i+r) <= rows || (j+c) <= cols );
  int aij;
  for (int ii=0; ii<r; ii++) {
    for (int jj=0; jj<c; jj++) {
      aij = this->getElement(i+ii, j+jj);
      this->setElement(i+ii,j+jj, aij+Sub.getElement(ii,jj));
    }
  }
  return *this;
}

//! Take sub-matrix from a given matrix
iMatrix iMatrix::getSubMatrix(int i, int j, int subRow, int subCol)
{
  assert (subRow > 0 && subCol > 0);
  iMatrix result(subRow, subCol);
  // check if submatrix belong to this matrix
  assert( (i+subRow) <= rows || (j+subCol) <= cols );
  int aij;
  for (int ii=0; ii<subRow; ii++) {
    for (int jj=0; jj<subCol; jj++) {
      aij = this->getElement(i+ii, j+jj);
      result.setElement(i+ii,j+jj, aij);
    }
  }
  return result;
}

//============================================= iMatrix operators
//!  Matrix B = Matrix self + Matrix A
iMatrix iMatrix::operator + (const iMatrix& A)
{
  assert(rows == A.getRows() && cols == A.getCols());

  iMatrix ans(rows,cols);
  ans = this->Add( A );
  return (ans);
}

//! Matrix B = Matrix self - Matrix A
iMatrix iMatrix::operator - (const iMatrix& A)
{
  assert(rows== A.getRows() && cols == A.getCols());

  iMatrix ans(rows,cols);
  ans = this->Sub( A );
  return (ans);
}

//--- Product of this matrix and an int
//! Matrix B = Matrix self * int s
iMatrix iMatrix::operator * (int s)
{
  iMatrix B(rows, cols);
  B = this->Dot( s);
  return ( B );
}

//--- Product of this Matrix and a Vector
//!   iVector b = self * iVector v
iVector iMatrix::operator * (const iVector& v) 
{
  // check compatibility of vectors
  assert(cols == v.getLen());
  iVector ans(rows);
  ans = this->Dot( v );
  return (ans);
}

//--- Product of this Matrix and Matrix A
//!   Matrix B = self * Matrix A
iMatrix iMatrix::operator * (const iMatrix& A)
{
  // check compatibility of matrices
  assert(cols == A.getRows());
  int Cols = A.getCols();
  iMatrix B( rows, Cols );  // result matrix: B = self * A
  B = this->Dot( A );
  return ( B );
}
//
//! Transpose of a SELF Matrix
//
iMatrix iMatrix::operator ~ ()
{
  iMatrix ans(cols, rows);              // ans will be transpose of SELF
  ans = this->Tra();
  // give solution
  return (ans);
}

//======================================= iMatrix operators - equivalent private functions
//!  Matrix B = Matrix self + Matrix A
iMatrix iMatrix::Add(const iMatrix& A)
{
  assert(rows == A.getRows() && cols == A.getCols());

  iMatrix ans(rows,cols);
  int Cij;
  for ( int i = 0; i < rows; i++ ) {
    for (int j = 0; j < cols; j++) {
      Cij = getElement(i,j) + A.getElement(i,j);
      ans.setElement(i,j,Cij);
    }
  }
  return (ans);
}

//! Matrix B = Matrix self - Matrix A
iMatrix iMatrix::Sub(const iMatrix& A)
{
  assert(rows== A.getRows() && cols == A.getCols());

  iMatrix ans(rows,cols);
  int Cij;
  for ( int i = 0; i < rows; i++ ) {
    for (int j = 0; j < cols; j++) {
      Cij = getElement(i,j) - A.getElement(i,j);
      ans.setElement(i,j,Cij);
    }
  }
  return (ans);
}

//--- Product of this matrix and an int
//! Matrix B = Matrix self * int s
iMatrix iMatrix::Dot(int s)
{
  iMatrix B(rows, cols);
  int a_ij;
  for ( int i = 0; i < rows; i++ ) {
    for ( int j = 0; j < cols; j++ ) {
      a_ij = getElement( i, j ) * s;
      B.setElement( i, j, a_ij );
    }
  }
  return ( B );
}

//--- Product of this Matrix and a Vector
//!   iVector b = self * iVector v
iVector iMatrix::Dot(const iVector& v)
{
  assert(cols == v.getLen());

  iVector ans(rows);
  int ai;
  for(int i=0; i<rows; i++) {
    ai = 0;
    for(int j=0; j<cols; j++) {
      ai += getElement(i,j) * v.getData(j);
    }
    ans.setData(i,ai);
  }
  return (ans);
}

//--- Product of this Matrix and Matrix A
//!   Matrix B = self * Matrix A
iMatrix iMatrix::Dot(const iMatrix& A)
{
  assert(cols == A.getRows());

  int Cols = A.getCols();
  iMatrix B(rows,Cols);            // resultant matrix: B = self * A
  int Bij;
  for ( int i = 0; i < rows; i++ ) {
    for ( int j = 0; j < Cols; j++ ) {
      Bij = 0;
      for ( int k = 0; k < cols; k++ ) {
        Bij += getElement( i, k ) * A.getElement( k, j );
      }
      B.setElement( i, j, Bij );
    }
  }
  return (B);
}

//============================== Transpose of a SELF Matrix
//
//! Transpose of a SELF Matrix
//
iMatrix iMatrix::Tra()
{
  int N = this->getRows();
  int M = this->getCols();             // SELF may be N x M MAtrix
  iMatrix ans(M,N);                    // ans will be transpose of SELF
  iVector* tmp = new iVector(M);       // row of SELF Matrix
  assert(tmp != NULL);
  for (int i = 0; i < N; i++) {
    *tmp = this->getRowData(i);
    ans.setColData(i,*tmp);
  }
  delete tmp;
  // give solution
  return (ans);
}

//----------------- console I/O
//
//! Print each row separately (with info about matrix)
void iMatrix::Print(void)
{
  cout << "iMatrix: Rows = " << rows << "Cols = " << cols << endl;
  for(int i=0; i<rows; i++) {
    m[i].rPrint();
  }
}

//! Print each row separately
void iMatrix::pPrint(void)
{
  for (int i=0; i<rows; i++) {
    for(int j=0; j<cols; j++)
      cout << m[i][j] << " ";
    cout << endl;
  }
}

//------------------ file I/O
//
//! read a Matrix from default file
iMatrix iMatrix::Read(void)
{
  string FileName = Name + ".txt";
  iMatrix ans(FileName.c_str());
  return (ans);
}

//! read a Matrix from a file
iMatrix iMatrix::ReadFrom(const char* file)
{
  iMatrix ans(file);
  return (ans);
}

//!  write Matrix to default file
void iMatrix::Write(void)
{
  int L = 20;                 // default number of entries in a line
  if (this->getCols() < L)
    L = this->getCols();
  string FileName = Name + ".txt";
  this->WriteTo(FileName.c_str(), L);
}

//!  write a matrix to a file
//
void iMatrix::WriteTo(const char* file, int L)
{
  int No_Lines;          // number of lines for ONE ROW
  int resto;             // left-over in the last line for the row
  div_t Res;
  Res.quot = 1;
  Res.rem = 0;
  if (cols <= L)
    L = cols;
  else if (cols > L)        // if number of cols is > number of entries in a row
    Res = div(cols, L);
  No_Lines = Res.quot;      // No_Lines with full number of mult entries
  resto    = Res.rem;       // ONE line with resto entries, all for just ONE ROW of data
  int J, data;
  string MatrixName = Name;
  //-- output to file
  out.open(file,ios::out);
  assert(out != NULL);
  //--- write the HEADER BLOCK
  out << "!================== iMatrix(" << rows << ", " << cols << ") ====================" << endl;
  out << "N= " << MatrixName << endl;
  out << "r= " << rows << "  c= " << cols << "  l= " << L  << endl;
  out << ":::" << endl;
  //--- write the DATA BLOCK
  int kmax;
  for (int i = 0; i < rows; i ++) {
    for (int j = 0; j <= No_Lines; j++) {
      kmax = L;
      if (j == No_Lines)
        kmax = resto;
      for (int k = 0; k < kmax; k++) {
        J = j * L + k;
        data = this->getElement(i,J);
        out << data << "   ";
      }
    }
    out << endl;
  }
  out.close();
}

//=============== Global iMatrix operator functions
//
//!  Multiply a matrix with a scalar: B = s * A
//
iMatrix operator * (int s, const iMatrix& A)
{
  iMatrix ans(A.getRows(),A.getCols());
  for(int i=0; i<ans.getRows(); i++) {
    ans.setRowData(i, s * A.getRowData(i));
  }
  return (ans);
}

//!  Multiply a matrix with a scalar: B = A * s
//
iMatrix operator * (const iMatrix& A, int s)
{
  iMatrix ans(A.getRows(),A.getCols());
  for(int i=0; i<ans.getRows(); i++) {
    ans.setRowData(i, s * A.getRowData(i));
  }
  return (ans);
}

//! Multiplication of a matrix and a vector: c = A * b
iVector operator * (const iMatrix& A, const iVector b)
{
  // check matrix - vector compatibility
  int M = A.getRows(), N = A.getCols();
  assert(N == b.getLen());
  iVector c(M);
  int A_ij, bj, ci;
  for (int i=0; i<M; i++) {
    ci = 0;
    for (int j=0; j<N; j++) {
      A_ij = A.getElement(i,j);
      bj = b.getData(j);
      ci += A_ij * bj;
    }
    c.setData(i,ci);
  }
  return (c);
}

//! Multiplication of two matrices: C = A * B
iMatrix operator * (const iMatrix& A, const iMatrix& B)
{
  // check matrix compatibility
  assert(A.getCols() == B.getRows());
  int M = A.getRows(), N = B.getCols();
  iMatrix C(M,N);
  int A_ik,B_kj,C_ij;              // samo C_ij treba
  // matrix product
  for (int i=0; i<M; i++) {
    for (int j=0; j<N; j++) {
      C_ij = 0;
      for (int k=0; k<A.getCols(); k++) {
        A_ik = A.getElement(i,k);
        B_kj = B.getElement(k,j);
        C_ij += A_ik * B_kj;
      }
      C.setElement(i,j,C_ij);
    }
  }
  return (C);
}

//! Addition of two matrices: C = A + B
iMatrix operator + (const iMatrix& A, const iMatrix& B)
{
  // check matrix compatibility
  int M = A.getRows(), N = A.getCols();
  assert(M == B.getRows() && N == B.getCols());
  iMatrix C(M,N);
  int A_ij,B_ij,C_ij;
  for (int i=0; i<M; i++) {
    for (int j=0; j<N; j++) {
      A_ij = A.getElement(i,j);
      B_ij = B.getElement(i,j);
      C_ij = A_ij + B_ij;
      C.setElement(i,j,C_ij);
    }
  }
  return (C);
}

//! Subtraction of two matrices: C = A - B
iMatrix operator - (const iMatrix& A, const iMatrix& B)
{
  // check matrix compatibility
  int M = A.getRows(), N = A.getCols();
  assert(M == B.getRows() && N == B.getCols());
  iMatrix C(M,N);
  int A_ij,B_ij,C_ij;
  for (int i=0; i<M; i++) {
    for (int j=0; j<N; j++) {
      A_ij = A.getElement(i,j);
      B_ij = B.getElement(i,j);
      C_ij = A_ij - B_ij;
      C.setElement(i,j,C_ij);
    }
  }
  return (C);
}

//! stream output
//
ostream& operator << (ostream& s, const iMatrix& M)
{
  for(int i=0; i<M.getRows(); i++) {
    s << M.getRowData(i);
  }
  return s;
}

}   // end of namespace smv_lib

//==================================================== end of iMATRIX CLASS
//====================================== end of file

